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I.  INTRODUCTION 

The  completed  research  herein  described  is  concerned  with 
theoretical /computational  modeling  of  mechanical  signal  transmission 
through  soil  media  from  subsurface  sources.  The  propagation  of  such 
signals  plays  a  central  role  in  detecting  or  sensing  dynamic  subsurface 
events.  A  source  below  the  surface  will  transmit  stress  waves  in  all 
directions.  These  waves  will  interact  with  the  surrounding  surface  and 
any  buried  objects,  and  will  produce  a  dynamic  wave  motion  field.  By 
incorporating  a  proper  measurement  field  array,  data  may  be  gathered  on 
the  field  to  predict  the  nature  of  the  source.  Mathematical  predictions 
of  both  surface  and  subsurface  wave  motion  patterns  would  act  as  a  guide 
for  optimal  sensor  arrays  and  logic  in  order  to  determine  these  source 
signatures. 

The  propagation  and  scattering  of  elastic  waves  by  cavities  and 
cracks  has  been  the  subject  of  numerous  investigations  during  the  past 
several  years.  Early  studies  have  dealt  with  the  cases  of  cavities  in 
infinite  media,  while  more  recent  work,  Sanchez-Sesma  [46],  Stone  et.  al. 
[50],  Mendelsohn  et.  al.  [34],  Achenbach  and  Brind  [3],  Datta  and 
El-Aklly  [18],  Shah  et.  al.  [47],  and  Abdul jabbar  et.  al .  [1],  has  been 
directed  at  scattering  problems  in  semi-infinite  domains.  For 
semi-infinite  regions  (near  field  soil  modeling)  the  presence  of  the  free 
surface  Influences  the  scattering  wave  field,  and  this  phenomena  becomes 
Important  in  a  wide  variety  of  dynamic  geomechanics  problems  including 
seismology,  blast  loading,  sensing,  mineral  exploration,  etc. 

Wave  propagation  problems  quickly  become  difficult  to  analyze  as  the 
boundary  geometry  becomes  complex.  At  this  point  one  usually  turns  to 
some  numerical  method  that  will  approximate  the  solution  for  this 
geometry.  The  methods  of  finite  differences  and  finite  elements  have 


become  very  popular  in  wave  propagation  studies.  These  methods  are 
sometimes  called  domain  methods  because  the  entire  field  of  the  physical 
problem  is  discretized  into  small  areas  or  cells.  For  example,  Reynolds 
[42]  and  Kelley  et.  al.  [25]  have  developed  finite  difference  codes  for 
two-dimensional  wave  propagation.  Other  finite  difference  and  finite 
element  work  can  be  found  In  Alder  [5]  which  reviews  a  broad  class  of 
problems.  When  modeling  wave  propagation,  these  domain  methods  have 
limitations  in  properly  representing  the  boundaries  of  the  body  under 
study  and  have  trouble  if  the  domains  are  large  or  infinite.  Bouchon  and 
Aki  in  a  series  of  articles  [8,9,10]  presented  an  alternate  numerical 
approach  called  the  discrete  wave  number  representation  method.  Ray 
methods  have  also  been  applied  for  approximate  answers  to  such  problems 
see  for  example  Achenbach,  Gautsen,  and  McMaken  [4], 

After  a  review  of  these  previous  techniques  one  finds  several 
difficulties  associated  with  each  method.  Consequently  a  somewhat  new 
approach  was  taken  in  which  a  numerical  technique  was  developed  based  on 
an  integral  formulation  of  the  elastodynamic  equations.  Integral 
formulations  in  elasticity  theory  have  a  long  history  dating  back  to  the 
19th  century.  However,  although  such  formulations  have  been  known  for 
quite  some  time,  their  application  for  developing  numerical  algorithms 
for  use  in  problem  solutions  is  relatively  new.  Recently  many 
investigators  have  been  employing  integral  formulations  with  numerical 
methods  (Including  finite  element  concepts)  to  solve  a  broad  class  of 
both  static  and  dynamic  elasticity  problems.  These  methods  have  become 
known  as  the  Boundary  Integral  Equation  Method  (BIE  or  BIEM)  or  the 
Boundary  Element  Method  (BEM).  The  recent  texts  by  Brebbia  [11],  Brebbia 
and  Walker  [12]  and  Banerjee  and  Butterfield  [6]  provide  general 
background  for  these  methods  as  they  are  applied  to  a  large  variety  of 


engineering  science  problems. 

With  regard  to  elastodynamic  problems,  Cruse  and  Rizzo  [17]  and  Cruse 
[16]  developed  a  BEM  approach  in  conjunction  with  Laplace 
transformation.  Niwa  et.  al.  [39,40]  and  Kobayashi  and  Nishimura  [27] 
used  BEM  methods  for  the  steady  state  solution  and  then  reconstructed  the 
transient  response  by  Fourier  synthesis.  Shaw  [48]  also  discusses  the 
steady  state  solution  using  a  Fourier  transform.  Cole,  Kossloff,  and 
Minster  [15],  Mansur  and  Brebbia  [32,33],  and  Misljenovic  [37,38]  have 
developed  a  time-dependent  formulation  which  eliminates  the  need  for 
Laplace  transformation  and/or  Fourier  synthesis.  Manolis  [31]  has 
recently  compared  these  three  BEM  methods  on  a  scattering  proble’  an 
infinite  medium,  and  found  that  all  three  methods  gave  good  resul' 

Using  the  BEM  with  transform  methods,  after  one  obtains  a  so.  •t-’  ,  in 
the  Laplace  or  Fourier  domain,  one  must  invert  it  to  solve  for  the  time 
domain  solution.  Transform  methods  also  require  the  solution  of  a  full 
matrix  and  do  not  take  advantage  of  the  causal  properties  that  the 
time-dependent  Green's  function  possesses.  The  direct  BEM  time  domain 
approach  has  economical  advantages  because  of  this,  as  it  permits  an 
explicit  stepping  scheme  which  will  result  in  a  banded  solution  matrix. 

Fundamental  to  the  application  of  the  Boundary  Element  Method  is  a 
Green's  function  (fundamental  solution)  to  the  governing  equations.  The 
ease  and  efficiency  of  the  method  is  connected  quite  closely  with  the 
nature  of  the  Green's  function  that  is  employed.  All  of  the  previously 
mentioned  studies  have  used  the  Green's  function  for  an  infinite  domain, 
thus  requiring  all  boundaries  to  be  discretized.  Kobayashi  and  Nishimura 
[28]  constructed  a  half-space  Green's  function  for  the  steady  state 
elastodynamic  case,  but  did  not  utilize  it.  It  has  been  demonstrated  by 
Telles  and  Brebbia  [51]  for  the  static  elasticity  case,  that  by  using  a 


half-space  Green's  function,  considerable  simplification  is  gained  for 
half-space  geomechanics  problems.  A  major  goal  of  this  research  was  to 
extend  Telles'  and  Brebbia's  idea  to  dynamic  wave  propagation  phenomena. 
This  has  involved  the  development  of  an  easily  computable  Green's 
function  for  a  half-space  domain,  and  its  imp! indentation  into  a  general 
purpose,  two-dimensional  boundary  element  numerical  code.  The  created 
numerical  methods  are  designed  for  half-space  problems  involving 
propagation  and  scattering  of  near  field  transient  elastic  waves. 
Solution  methods  for  both  anti-plane  strain  (out-of-plane  motions)  and 
plane  strain  (in-plane  motions)  are  developed. 


II.  BASIC  ELASTODYNAMIC  THEORY 
2.1  Basic  Equations  of  Elastodynamics 

The  basic  variables  In  dynamic  elasticity  theory  are  the  stresses 
ojj,  the  strains  e^-  and  the  displacements  u7-.  The  general  field 
equations  relate  these  variables  to  one  another  and  provide  the  basis  of 
the  theory  (see  Achenbach  [2],  Graff  [23]  or  Ewing  et.al.  [21]).  From 
kinematics  the  strain-displacement  relations  are  given  by 


e. .  =  7T  ( u .  .  +  u.  •) 

i  j  2  v  i  ,j  j  ,v 


(2.1) 


Due  to  the  conservation  of  linear  momentum  we  can  express  Newton's  law  as 
the  equation  of  motion  for  a  conti num, 


U.J 


pfi  - 


pu. 


(2.2) 


where  f^  is  the  body  force  density  and  p  is  the  mass  density.  The 
constitutive  law  for  a  linear  elastic  isotropic  material  is  given  by 


Hooke's  law  as 


°1J  =  Xekksij  +  2Meij’  ,2-3) 

where  x  and  u  are  elastic  constants. 

Combining  equations  (2.1),  (2.2),  (2.3),  we  can  eliminate  the 
stresses  and  strains  to  get  Navier's  equation  of  motion  in  terms  of 
displ acement 


pui  ,kk  +  (A+p)uk,ki  +  pfi  =  p“i  (2*4) 

There  are  normally  three  types  of  boundary  conditions  prescribed  on 
the  boundary  surface  of  a  body. 

1)  A  distribution  of  surface  forces  or  tractions  _t,  on  the  surface 
of  the  body,  where 

‘i  =  °ij"i 

with  n^  being  the  unit  normal  vector. 

2)  A  distribution  of  displacements  on  the  surface  of  the  body. 

3)  A  distribution  of  surface  forces  on  a  portion  of  the  boundary 
surface  and  displacements  on  the  remaining  portion. 

These  along  with  initial  conditions  on  the  displacement  and  velocity 
provide  the  necessary  conditions  which  when  coupled  with  the  governing 
equations  produce  a  unique  solution.  The  proof  of  this  can  be  found  in 
Wheeler  and  Sternberg  [52]  or  Sokol nikoff  [49], 


2.2  Dilatation  and  Shear  Wave  Motion 


If  we  take  the  divergence  of  the  equation  of  motion  (2.4)  without 
body  forces,  we  get  an  equation  associated  with  irrotational  wave  motion, 
i.e. 


u 


1  " 

i  ,ikk  “  ~1  ui  ,i  ’ 
C1 


(2.5) 


where  u^  .  is  the  dilatation  (volume  changing  motion)  and  c-j  = 

[(  X  +  2  y  )/  This  irrotational  or  dilatational  wave  motion  is 

associated  with  what  is  commonly  called  a  P-wave  which  propagates  with 
velocity  c-,. 

Taking  the  curl  of  equation  (2.4)  without  body  forces,  we  get  an 
equation  associted  with  rotational  wave  motion. 


u 


i  ,kk 


(2.6) 


where  c2  ■  (y/p  )^.  This  rotational  (volume  preserving  motion)  or 
shear  vave  motion  Is  commonly  refered  to  as  an  S-wave  and  propagates  with 
velocity  c2. 

In  an  infinite  medium  an  Internal  input  disturbance  will  in  general 
produce  both  dilatation  and  shear  motions.  Thus  a  general  disturbance 
will  propagate  to  remote  points  as  two  different  wave  motions. 


2.3  Helmholtz  Decomposition  Theorem 

The  Helmholtz  Decomposition  Theorem  states  that  a  vector  field  u  that 
is  sufficiently  smooth  i.e.  satisfies  certain  continuity  requirements, 
may  be  written  as  the  sum  of  the  gradient  of  a  scalar  potential  function 
plus  the  curl  of  a  vector  function,  i.e. 


U  =  V<t>+Vx^, 


(2.7) 


with  the  condition 

V*£  =  0  . 

Using  this  form  of  £  as  a  displacement  vector  and  substituting  into 
equation  (2.4)  without  body  forces,  yields 

V^<p  =  — — <p  »  V^lp  =  — — jr  ip  .  (2.8) 

C1  c2 

Note  that  the  displacement  potential  <p  Is  associated  with  the  dilatation 
wave  and  ^  is  associated  with  the  shear  wave. 

This  decomposition  Is  useful  for  many  elastodynamic  problems,  since 
the  complex  equation  of  motion  (2.4)  is  reduced  to  two  simple  uncoupled 
wave  equations.  It  should  be  mentioned  here  however,  that  the  two 
potentials  are  coupled  through  the  boundary  conditions  since  the 
displacements  and  stresses  depend  on  both  $  and  V  . 

2.4  Two  Dimensional  Elastodynamlcs 
Plane  Strain 

For  bodies  that  are  very  long  in  one  direction  and  loaded  in  a  plane 
perpendicular  to  this  direction  we  can  make  certain  approximations  that 
will  simplify  the  general  equation  of  motion  (2.4).  This  case  is  known 
as  plane  strain  and  the  displacements  are  of  the  form 


u  *  “(X^xg.t) 

v*v(Xl,x2,t)  (2.9) 

w  ■  0 


Note  that  this  corresponds  to  in-plane  motions. 


-8- 


The  non-zero  strain  components  are  given  by 

3u 


ex  “  3x 


1 


3v_ 

3X, 


1_  /  3u  ,  3v  \ 
exy  ‘  2  ' 3x£  ' 

and  the  non-zero  stress  components  are 

°x  =  ^ex+ey)+2uex 
°y  =  X(ex+ey)+2yey 
az  =  X(ex+ey) 

Txy  =  2yexy- 


(2.10) 


(2.11) 


The  equation  of  motion  in  terms  of  displacement  becomes 


u  v^u+  ( A+U ) )+p£  =  P]j 


(2.12) 


where 


Waves  associated  with  this  equation  will  either  be  dilatation  waves 
(P-waves)  or  vertically  polarized  shear  waves  (SV-waves). 

Plane  Stress 

For  bodies  that  are  very  thin  in  one  direction  and  loaded  in  a  plane 
perpendicular  to  this  direction  we  can  again  make  certain  approximations 
that  will  simplify  the  basic  equation  of  motion.  This  case  is  refered  to 
as  plane  stress  and  the  stress  components  are  taken  to  be 


-9- 


°x  =  °x^xl ,x2*t^ 

°y  =  ay^X1 »x2‘^) 

Txy  =  Txy(xl ,x2,t} 
az  =  Txz  =  Tyz  =  0 


The  non-zero  strain  components  will  then  be 

ex  "  <Vwy)/E 

ey  ■  (Oy+VOyl/E 

e  =  (o  +a  )/E 
z  x  y" 

e  =  t  /2u, 
xy  xy 


(2.13) 


(2.14) 


while  the  equation  of  motion  In  terms  of  displacement  reduces  to 


yV2u  +  2 ( Z(V*u)+pf  =  pu 


(2.15) 


where  E  and  v  are  elastic  constants. 


Anti -Plane  Strain 

Bodies  that  are  very  long  in  one  direction  and  loaded  parallel  to 
that  direction  can  again  be  modeled  by  making  certain  approximations  that 
will  simplify  the  basic  three  dimensional  equations  of  motion.  This  case 
is  normally  referred  to  as  anti-plane  strain  and  the  displacements  are  of 
the  form 


u  *  0 
v  ■  0 

w  *  w(x1 ,x2,t). 

This  case  corresponds  to  out-of-plane  motions. 


(2.16) 


The  non- zero  strain  components  follow  from  (2.1)  to  be 


eyz  =  3x2  (2.17) 

-  9w 
exz  "  3x1 

and  the  non-zero  stress  components,  from  Hook's  law,  are 


a  =  2ue 
zx  ^  zx 


(2.18) 


a  =  2ue 
zy  M  zy 


For  this  case  the  equation  of  motion  in  terms  of  displacement  reduces 


to 


uv2w+pf  =  pw  (2.19) 

Waves  associated  with  this  anti-plane  case  will  be  horizontally  polarized 
shear  waves  (SH-waves). 

2.5  Ray  Theory 

The  purpose  of  this  section  is  to  review  how  plane  P-,SV-,  and 
SH-waves  reflect  from  stress  free  plane  surfaces.  The  geometrical  nature 
of  these  incident  and  reflecting  waves  is  sometimes  called  ray  theory. 
For  plane  strain,  the  free  surface  will  produce  a  coupling  (mode 
conversion)  between  the  SV-wave  and  the  P-wave  i.e.,  the  incident  P-wave 
or  SV-wave  will  produce  both  reflected  P-waves  and  SV- waves. 

All  cases  are  related  through  Figure  2.1  which  illustrates  the  basic 
ray  patterns.  The  incident  wave  (either  P  or  S)  has  an  incident  angle  Q0» 
measured  from  the  surface  normal.  The  reflected  P-wave  has  ancle  ,9-| , 
while  the  reflected  S-wave  has  angle  Og. 


Fig.  1  Ray  Theory 


For  the  anti -plane  strain  case,  there  will  be  no  mode  conversion  and 
hence  only  a  single  reflected  SH-wave  will  be  present  with  the  incident 
and  reflected  angles  being  equal,  i.e. 


For  plane  strain  we  have  two  possible  types  of  incident  waves  the 
P-wave  and  the  SV-wave.  In  both  cases  there  is  in  general  a  mode 
conversion  where  either  the  P-wave  or  SV-wave  will  upon  reflection  become 
part  P-wave  and  part  SV-wave;  see  Figure  2.1. 

For  an  incident  P-wave, 


9 


1 


(2.21) 


c2$in  0Q  =  c^sinSg, 


(2.22) 


while  for  an  incident  SV-wave 


e 


2 


(2.23) 


CgSine-i  *  c-]Sin90 


(2.24) 


Note  equation  (2.23)  is  real-valued  only  if  e<sin_1  (c2/c-| )  otherwise  there 
is  a  mode  conversion  of  the  SV-wave  to  a  surface  wave. 

The  previous  equations  provide  information  necessary  to  compute 
directly  the  angles  and  ray  path  distances  of  reflected  and  incident 
plane  waves,  except  in  the  case  of  the  reflected  SV-wave  from  an  incident 
P-wave  and  the  reflected  P-wave  from  an  incident  SV-wave.  To  find  the 
ray  path  distances  for  these  two  exceptions  we  must  turn  to  numerical 
techniques. 

Figure  2.2  illustrates  a  typical  source/ receiver  pair  in  plane 
strain,  where  the  source  is  at  (0,h)  and  the  receiver  is  at  (x,y). 
Consider  a  ray  path  SPR  as  shown,  depicting  an  incident  P-wave  traveling 
a  distance  c,  and  a  reflected  SV-wave  traveling  a  distance  d.  With  x  = 
a+b,  these  ray  paths  c  and  d  can  be  found  from  two  simultaneous  nonlinear 
equations 


2  r  Kxd  -,2  _  2 
h  +[x  •  c+Kcr  '  C 


2  -Kxd  n2  .2 
*  +C^«d]  =  d  • 


(2.25) 


where,  K  =  c2/c-|.  These  equations  can  be  solved  by  employing  the 
Newton  method  see  Rice  [45],  Ketter  and  Sherwood  [26]  and  Berezin  and 
Zidkov  [7]. 
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Fig.  2.2  Geometry  of  Incident  P-Wave 
and  Reflected  SV-Wave 


Figure  2.3  illustrates  another  typical  source/receiver  pair  in  plane 
strain,  where  the  SV-source  Is  at  (0,h)  and  the  receiver  Is  at  (x,y). 
Consider  now  a  ray  path  SPR  as  shown,  where  an  Incident  SV-wave  travels  a 
distance  c,  and  a  reflected  P-wave  travels  a  distance  d.  With  x  =  a+b, 
these  ray  paths  can  be  found  from  two  simultaneous  equations 


h^+r  _xd-j2  2 
n  L*  Kc+dJ  c 

2^r  xd  n2  .2 


(2.26) 


These  equations  can  again  be  solved  by  employing  the  Newton  method. 

Figure  2.4  Illustrates  the  case  of  waves  with  equal  angles  of 
Incidence  and  reflection.  This  occurs  for  the  SH,  P-P  and  SV-SV  cases. 
Again  the  source  Is  at  (0,h)  and  the  receiver  is  at  (x,y),  and  c  and  d 
are  the  ray  paths  of  incidence  and  reflection.  For  this  case 
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Fig.  2.3  Geometry  of  Incident  SV-Wave 
and  Reflected  P-Wave 


c  = 


x  2  1/2 

h[l  +  (t5^)2] 


(2.27) 


2.6  Integral  Representations 

For  the  general  elastodynamic  problem,  we  may  obtain  an  integral 
equation  from  the  dynamic  reciprocal  theorem  which  is  the  direct 
extension  of  Betti's  reciprocal  theorem  in  elastostatics.  The  theorem 
states  that  two  unrelated  elastodynamic  states  consisting  of  body  forces, 
surface  tractions  and  displacements,  initial  conditions  and  velocities, 
such  as  (fj ,t. ,u.,u^,v!}  and (f^ ,t. ,u. ,u^,v! )  which  exist  on  the  same  body  S 


9 


RaflMtM 

SH-V«r* 


Incident 


Fig.  2.4  Geometry  of  Incident  SH-Wave 
and  Reflected  SH-Wave 


bounded  by  a  surface  3S  for  t>  C,  are  related  by 

/  C(t1*u1.)(x,t)]dS  +  /pHfj^UjMx.tJ+Vj * (x)u^(x*t) 
3S  S 

3u. 

+  ^'(x)  (x,t)]dv  «  /[ti*u1)(x,t)]dS 

3S 

3u.(x,t) 

+  ;c[(fi*ui)(x,t)  +  v!(x}ui(x,t)+ui ' (x)  —  ]dv 


where  *  means  the  convolution  integral  defined  by 


if*9)(x,t) 


f  f(x,t-x)g{x,t)dT  ,  t  >  0 
o 

0  ,  t  <  0  . 


If  we  choose  one  of  the  elastodynamic  states  say  (f^ .t^ .u^ ,u!  ^ 
be  that  of  the  physical  problem  and  let  the  other  state  (f1»t.j,ui,u 


that  generated  by  a  unit  concentrated  body  force  In  the  same  body,  we 
obtain  the  Integral  expression 


00 

cu(ro.t0)=  |{Si(r,t)-u(r,t)j£#n}dS(E)dt, 
0  35 


(2.30) 


where  we  have  neglected  the  body  force  f1  and  taken  zero  Initial 
conditions.  The  coefficient  c  Is  given  by 


c 


1/2,  r  e  smooth  boundary  3S 
1,  r  e  Interior  of  as 
0,  r  e  exterior  of  3S 


(2.31) 


d* .  _ 


G^G(r,t;r0,t0)  is  the  displacement  Green's  function,  K  Is  the  stress 
field  associated  with  G,  x  and  u  are  the  surface  tractions  and 
displacements  on  the  boundary  3S  ,  Io  and  ^  are  the  space  and  time 
positions  of  the  source  (associated  with  the  Green's  function),  r  and  t 
are  the  space  and  time  positions  of  the  receiver,  and  Is  the  unit 
normal  vector  to  the  surface.  Subscripts  "o"  refer  to  source 
quantities.  Details  on  the  Green's  functions  will  be  provided  later  In 
the  report. 

For  the  ant1>plane  strain  case,  the  Integral  relation  may  be  written 
as 


cw(r,t) 


3 }  j  {^(ro’V-^vV^dstroJdV 

0  dS 


(2.32) 


where  we  now  have  a  scalar  equation. 


Ill  THE  BOUNDARY  ELEMENT  METHOD 


3.1  Principal  Advantages 

Recently  a  new  numerical  method  has  been  developed  called  the 
Boundary  Integral  Equation  Method.  The  basic  method  when  combined  with 
finite  element  discretization  concepts  Is  commonly  refered  to  as  the 
Boundary  Element  Method.  The  method  Is  applicable  to  a  wide  variety  of 
field  problems  In  engineering  science  and  has  some  advantages  over  the 
finite  element  and  finite  difference  methods  In  several  situations. 

The  principal  advantages  of  the  BEM  method  lie  primarily  In  the 
reduction  of  the  number  of  spatial  dimensions  of  the  problem  by  one,  and 
by  formulating  the  problem  directly  In  terms  of  the  boundary  values. 
Since  discretization  Is  carried  out  only  on  the  boundary  of  the  body 
under  study,  the  number  of  unknowns  In  the  numerical  problem  Is  reduced, 
and  the  need  of  creating  space-filling,  three-dimensional  grids  Is 

eliminated.  For  many  cases,  problems  Involving  Infinite  domains  can  be 
handled  In  a  very  direct  and  simple  manner  using  appropriate  Green's 
functions.  For  Infinite  or  semi-infinite  problems  there  will  be 
considerably  fewer  elements  and  thus  a  smaller  system  of  equations  need 
be  solved. 

Another  advantage  over  domain  methods  has  to  do  with  the  solution 
variables  Inside  the  body.  These  variables  In  domain  methods  will  often 
show  unrealistic  jumps  In  values  between  nodal  points  and  will  not  vary 
continuously.  In  the  BEM  method  these  variables  will  vary  continuously 
In  the  body  as  the  only  approximations  of  geometry  occur  on  the 

boundary.  In  many  engineering  problems  the  Important  part  of  the 

solution  will  be  on  the  boundary  or  at  a  few  selective  points  In  the 

Interior,  with  domain  methods  the  solution  Is  found  at  all  nodal  values, 
while  In  the  BEM  method  the  solution  Is  found  for  boundary  nodes  and  only 


those  Interior  points  selected  by  the  user.  As  In  the  finite  difference 
and  finite  element  methods,  the  problem  Is  reduced  to  a  system  of 
simultaneous  equations  to  be  solved.  While  In  the  domain  methods  the 
system  coefficient  matrix  Is  usually  banded,  using  the  BEM  this  matrix 
can  be  fully  populated. 

3.2  Method  Development 

The  development  of  the  BEM  method  for  a  general  problem  Is  described 
in  the  flow  chart  In  Figure  3.1.  The  method  Is  started  with  the  use  of 


Fig.  3.1  Flow  Chart  of  Boundary  Element  Method 


potential  theory  or  weighted  residual  methods  which  when  combined  with  a 
Green's  function,  form  the  boundary  Integral  equation.  The  use  of  an 
efficient  Green's  funclton  Is  central  to  the  method.  There  are  multiple 
choices  for  a  Green's  function,  all  of  which  will  satisfy  the  governing 
equation  for  the  problem,  but  will  vary  according  to  the  boundary 
conditions  on  the  body.  Again  referlng  to  Figure  3.1,  the  boundary 


Integral  equation  Is  then  discretized  using  established  finite  element 
techniques  to  form  a  system  of  simultaneous  linear  algebraic  equations. 
The  solution  to  these  equations  will  yield  the  unknown  boundary  values  of 
the  problem.  These  boundary  values  can  In  turn  be  used  to  find  the 
solution  at  any  Interior  points. 

The  basis  for  an  integral  equation  formulation  In  elastodynamics  was 
first  developed  by  Kupradze  [29],  de  Hoop  [19]  and  Wheeler  and  Sternberg 
[52],  while  Erlngen  and  Suhubl  [20]  provide  a  concise  review  of  this 
classical  work.  The  Integral  equation  formulation  Is  the  dynamic 
extension  of  Betti's  reciprocal  theorem  which  relates  two  different 
dynamic  states  for  a  given  elastic  body,  see  equation  (2.28).  An 
apparent  choice  for  one  of  the  two  elastodynamlc  states  Is  the  desired 
solution  to  the  problem,  while  the  other  state  Is  normally  chosen  to  be 
an  appropriate  Green's  function  that  satisfies  the  governing  equations  of 
motion  for  the  problem  under  study.  Using  these  choices,  the  dynamic 
equivalent  to  Somlgllana's  Identity  can  be  derived,  e.g.  equation  (2.30), 
and  through  a  limiting  process  a  boundary  Integral  equation  Is  developed 
which  forms  the  basis  of  the  BEM  method.  This  process  was  outlined  In 
section  2.6,  and  developes  what  Is  called  a  direct  BEM  approach,  see 
Mendel  son  [35]. 

3.3  Half-Space  Green's  Function  for  Geomechanics  Problems 

If  one  were  to  model  a  typical  geomechanics  problem  of  a  burled  hole 
using  a  finite  element  technique,  a  typical  mesh  might  look  like  Figure 
3.2.  As  can  be  seen  In  this  figure  there  are  problems  in  determining  how 
far  out  to  run  the  mesh  and  what  boundary  conditions  to  use  on  the  outer 
elements.  The  BEM  method  will  eliminate  the  need  for  a  space  filling 


Fig.  3.2  Finite  Element  Interior  Discretization 
for  a  Buried  Hole 


mesh  as  can  be  seen  in  Figure  3.3a  where  this  discretization  is  a  result 
of  using  the  BEM  with  an  infinite  space  Green's  function.  Again  however, 
problems  appear  as  to  how  far  out  to  run  the  boundary  mesh  and  what 


a.  Infinite  Space  Green's  Function  b.  Half-Space  Green's  Function 

Fig.  3.3  BEM  Discretization  for  Buried  Hole 


boundary  conditions  to  use  on  the  outer  elements.  These  problems  can  be 
eliminated  with  the  use  of  a  half-space  Green's  function  as  shown  In 
Figure  3.3b,  where  now  there  Is  no  need  to  discretize  the  free  surface 
since  that  boundary  condition  is  automatically  satisfied  by  the 
half-space  Green's  function.  This  example  brings  out  the  importance  of 
the  choice  of  the  Green's  function,  and  illustrates  the  usefulness  of  a 


half-space  form  for  geomechanics  problems.  The  employment  of  such  a 
half-space  Green's  function  was  a  major  part  of  this  research  program. 

TV  GREEN'S  FUNCTIONS 

4.1  Introduction 

The  Green's  function  Is  the  singular  solution  of  a  differential 
equation  subjected  to  an  Impulsive  source.  In  the  case  of  elastodynamlcs 
the  governing  equation  Is  Navler's  equation  of  motion  (2.4),  and  the 
source  Is  taken  as  a  concentrated  Impulsive  body  force  loading.  The 
Green's  function  may  be  developed  for  regions  of  Infinite  extent  or  for 
regions  of  finite  size  with  particular  boundary  conditions.  As 
demonstrated  In  section  3.3,  the  development  of  a  Green’s  function  for  a 
semi-infinite  domain  Is  most  desired  for  geomechanics  applications. 

Elastodynamlc  problems  for  half-spaces  subject  to  Impulsive  burled 
point  loads  or  other  sources  have  been  studied  for  quite  some  time,  see 
e.g.  Lamb  [30],  Hudson  [24],  Erlngen  and  Suhubl  [20],  Payton  [41],  and 
Garvin  [22].  These  types  of  problems  are  directly  related  to  fundamental 
singular  solutions  (l.e.  Green's  functions)  for  the  half-space  geometry. 
For  the  anti-plane  strain  case,  such  solutions  can  be  found  In  Cole 
et.al.  [15],  Achenbach  [2],  Rice  and  Sadd  [44],  and  Hudson  [24], 
Singular  force  solutions  for  the  plane  strain  case  are  complex  and  are 
obtained  by  integral  transform  methods  with  the  inversion  being 
accomplished  by  employing  some  version  of  the  Cagnlard  [13]  technique. 
Because  of  mathematical  difficulties,  most  of  these  solutions  are 
formally  carried  out  only  for  surface  points  or  other  special  points  in 
the  medium. 
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Wlth  a  fundamental  Green's  function  for  a  half-space  geometry,  the 
BIE  method  can  be  used  quite  effectively  to  handle  many  dynamic 
geomechanics  problems,  see  for  example  Rice  and  Sadd  [43,44].  However  In 
order  to  use  a  fundamental  solution  In  this  way.  It  must  be  computable  at 
all  points  In  the  region.  Hence  the  development  of  easily  computable 
Green's  functions  Is  very  Important  In  constructing  efficient  numerical 
BEM  codes  for  wave  propagation. 

We  present  here  the  development  of  an  elastodynamlc  Green's  function 
for  a  two-dimensional.  Isotropic  half-space.  For  anti-plane  strain,  we 
consider  the  burled  singular  impulsive  line  load  and  construct  a 
half-space  Green's  function  from  the  Infinite  space  Green's  function  by 
the  method  of  Images.  For  plane  strain,  we  address  the  burled  singular 
impulsive  line  load  as  formulated  by  Hudson  [24].  Starting  with  the 
formal  solution,  a  numerical  root  finding  scheme  Is  employed  to  compute 
terms  which  normally  present  problems  for  points  below  the  free  surface. 

4.2  Anti -Plane  Strain  Green's  Function 

Rewriting  the  anti-plane  strain  equation  of  motion  equation  (2.19) 

yV2w+pf  =  pw,  (4.1) 

if  we  let  pf  be  a  concentrated  impulsive  out  of  plane  line  load  at 
rQ  and  tQ,  w  will  now  become  the  anti-plane  strain  Green's  function 
G.  For  this  case  equation  (4.1)  becomes 

yV2G-pG  «  -<S(t-tQ)6(r-r0).  (4.2) 


For  the  Infinite  space  domain,  G  Is  given  by 
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1  H(t-t  -R/Cp) 

s(r-‘ir0.‘0>  •  ss  [(t.t#)4.R2/c24]1/4 


(4.3) 


where  R  *  | .  Hand  H  is  the  Heaviside  function. 

The  stress  field  associated  with  this  Green's  function  is  computed 
through  Hooke's  law  and  is  given  by 


D  H(t-t  -R/Cp)  »sp 

K(r,t;r  ,t)  =  P  P  5  p  3/2  "sn  * 

“  "°  2irc22  C(t-t0)2-R2/c22]3/2 


(4.4) 


For  the  anti-plane  strain  half-space  problem,  the  method  of  Images 
can  be  used  to  construct  a  half-space  Green's  function  from  the  infinite 
space  Green's  function  through  the  superposition  of  a  virtual  Infinite 
space  function  of  equal  amplitude  with  an  image  source  point  relative  to 
the  free  surface  as  shown  In  Figure  4.1.  The  associated  stress  field  can 

IMAGE  SOURCE  NODE 
(*i.-*a) 


RECEIVER  NODE 


(*i'*a) 

SOURCE  NODE 


Fig.  4.1  Geometry  of  Source  and  Image  for  the  Construction 
of  a  Half-Space  Green's  Function 


be  derived  as  in  (4.4).  The  results  are 


r-r+G-  1  f  H(t‘VRo/c2) 

G*60+6i 


H(t-t0-R1/c2) 


C(t-t0)2-R02/c2231/2  C(t-t0)2-RT2/c22:1/2 


}  (4.5) 


K=K +K,  = 
— -o  —I 


Vl'-W'z1 

~7  1 - jn — 5 - 7 

2^2  C  (t-tQ )  “Rq  /c2^ 


•{ 


3R_  RiH(t-t  -Ri/c«)  9Ri 
~~T77  TfT-  + - 7 - 2 - 2~ 577  Tn~ ^  ( 4  •  ® ) 

r72  3-  c ( t-t0 j^r, 2/c22]3/z 


where  the  position  vectors  RQ  and  R1  are  defined  in  Figure  4.1.  This 
simple  superposition  produces  a  stress  free  surface  at  x2  =  0  as 

required.  The  advantage  of  using  the  half-space  form  in  the  BEM  method 
is  that  the  boundary  of  the  free  surface  will  not  have  to  be  discretized, 
thus  simplifying  the  problem. 

4.3  Plane  Strain  Green's  Function 

The  plane  strain  equations  of  motion  (2.12)  in  indicia!  notation  are 

uukin  ♦  (Xtu)Uf>1k  +  pfk  =  suk  (4.7) 

If  we  let  pfk  be  a  concentrated  in-plane  line  load  located  at  =  0, 

x2=h  and  u^G^ej  where  ej  is  a  unit  vector  in  the  x-j  ,x2 

plane,  then  G^j  is  interpreted  as  the  two-dimensional  Green's 

function.  This  Green's  function  G^j  evaluates  the  displacement  field 
Uj  due  to  a  concentrated  unit  force  in  the  e^  direction.  Equation 

(4.7)  then  becomes 

yGij,kk  +  (X+y)Gkj,ik  '  pGij  =  6(t)6(x-|  )6(x2-h)6...  (4.8) 
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Introducing  displacement  potentials  and  using  Fourier  transformation 
methods  with  the  Cagnlard  [13]  Inversion  technique,  Hudson  [24]  develops 
a  solution  to  (4.8)  for  a  stress  free  surface  with  zero  Initial 
conditions.  The  reflected  portion  of  the  wave  field  Is  found  to  be 


6?j  *  2iS  +  H(t-tsp)ufP  *  H(t-tps>u£ 


+  H(t-r"/c2)U^  +  H(t-tSPS)H(r7c2-t)U^S], 


(4.9) 


where  t  Is  the  actual  time,  r“  Is  the  distance  between  the  image  source 
and  receiver,  and  r'  Is  the  distance  between  the  source  and  receiver. 


PP 

U1d  is  the  displacement  tensor  for  an  Incident  P-wave  producing  a 


reflected  P-wave, 


s\  -is 


S(s)  ds  i 

R(sT  dt '  ’ 


(4.10) 


Is  the  displacement  tensor  for  an  Incident  P-wave  producing  a 
reflected  S-wave, 


,sns  Vs 

-  Re  { 

1 J  O 

-s  isru 


Q(s)  ds  i 
R(sT  dt ;  ’ 


(4.11) 


Ujj  is  the  displacement  tensor  for  an  Incident  S-wave  producing  a 
reflected  S-wave, 


Ui j  “  Re  { 


n«  -is 


is  s  /n 


S(s)  ds  i 

RfiT  dt  *  ' 


(4.12) 


SP 

is  the  displacement  tensor  for  an  Incident  S-wave  yielding  a 


reflected  P- wave. 
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SPS 

U.jj  Is  the  displacement  tensor  for  a  Incident  S-wave,  yielding  a 
surface  P-wave,  and  a  reflected  S-wave, 
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(4.14) 


where 


•  J 


na  •  (s2-!*,2)1'2 
n6  -  (s2-1/c22)'/2 

Q(s)  *  4is(2s2-l/c22) 

(4.15) 

L- 

R(s)  =  (2s2-l/c22)2-4s2nang 

S(s)  =  4is(2s2-l/c22) 

From  the  Cagnlard  method,  there  Is  a  functional 

relationship  between  s 

and  t  for  each  contribution,  l.e. 

. 

For  the  PP  term:  t  *  sx^j  +  ina(h+x2) 

(4.16) 

• 

-  :•  -t-i 

For  the  PS  term:  t  «  sx1  +  i(ngX2+nah) 

(4.17) 

For  the  SS  term:  t  =  sx1  +  ing(h+x2) 

(4.18) 

'  ‘  ; 

For  the  SP  term:  t  *  sx,  +  i(n  x,+nQh) 

1  CL  C  p 

(4.19) 

ft _ _ 

Equation  (4.9)  corresponds  to  Hudson's  equations  (7.66)  and  (7.67)  where 

we  have  used  a  more  compact  notation.  The  five  reflection  terms  In  (4.9)  1 

correspond  to  the  five  standard  reflected  wave  types  as  shown  In  Figure 


4.2. 


receiver 


Fig.  4.2  The  Ray  Paths  from  Source  to  Receiver 

As  is  common  in  problems  of  this  type,  some  of  the  terms  in  relation 

(4.9)  cannot  be  evaluated  In  closed  form  for  a  general  field  point  below 

SP  PS 

the  boundary  surface.  For  this  case,  the  terms  and  are  the 

ones  which  posses  this  problem.  Through  a  Cagnlard  transformation,  the 

quantity  s  must  be  determined  by  Inverting  the  relations  In  equations 

(4.16)-(4.19).  Equations  (4.17)  and  (4.19)  unfortunately  do  not  lend 

themselves  to  simple  Inversions,  and  consequently  a  Newton  root  finding 

method  for  analytic  functions  was  employed,  see  Rice  [45]  for  details. 

For  specified  values  of  x^t  Xg  and  t,  the  numerical  routine  computes 

s  such  that  these  equations  will  be  satisfied.  The  derivative  terms  in 

equations  (4.11)  and  (4.13)  are  computed  by  first  employing  Implicit 

differentiation  of  relations  (4.17)  and  (4.19)  and  then  using  the  root 

SP 

finding  values.  This  then  allows  direct  computation  of  the  and 
PS 

Ujj  terms. 

The  stress  field  associated  with  this  Green's  function  is  computed 


through  Hooke's  law  to  be 


K1jk  ■  [«,jVs+v(«irV5isV)]Grk.s 


(4.20) 


4.4  Results  for  the  Plane  Strain  Case 

Specific  results  for  the  plane  strain  case  were  Investigated  for  five 
points  within  the  half-space  as  shown  In  Figures  4. 3-4. 6.  The  coordinate 
points  are  (0,5),  (3,5),  (5,5),  (5,3),  (5,0)  all  In  cm.  The  Figures  plot 
the  non-dimensional  Green's  function  2c2G.j ^  rcpR  versus  dimensionless 
time  c2t/R,  where  R  Is  the  distance  from  the  source  to  the  receiver. 
The  material  properties  chosen  for  use  correspond  to  c-|  *  5800  m/s, 
c2  -  3350  m/s  and  p  »  2.5xl03  kg/m3,  and  h  »  10cm. 

For  the  point  (0,5)  In  Figures  4.3  and  4.4  the  solid  line  is  to  be 
compared  with  the  dots  which  represent  the  analytical  solution  by  Payton 
[41].  Also  for  the  point  (5,0)  In  Figures  4.5  and  4.6,  the  solid  line  Is 
compared  with  the  dots  which  are  the  analytical  solution  results  by 
Erlngen  and  Suhubl  [20].  As  can  be  seen  our  results  compare  quite  well 
to  these  analytical  solutions  for  the  special  points  considered. 

The  surface  results  at  (5,0)  as  shown  In  Figures  4.5  and  4.6 
Initially  presented  a  problem  by  not  decaying  properly  after  all  the  wave 
fronts  had  passed.  This  problem  did  not  occur  for  the  subsurface  output 
in  Figures  4.3  and  4.4.  This  decay  problem  was  a  result  of  a  numerical 
Instability  that  occurred  when  the  Individual  reflection  components  were 
very  large  compared  to  the  total  sum.  Each  total  displacement  component 
should  undergo  monotonic  decay  with  Increasing  time  after  all  the  wave 
fronts  have  gone  by. 


In  order  to  resolve  this  Instability  problem  a  smooting  operation  was 
employed.  After  all  the  wave  fronts  have  arrived.  If  Instability  (I.e. 
growth)  was  detected,  the  response  was  forced  to  decay  at  the  same  rate 
as  an  Infinite  space  source  problem  with  a  similar  receiver  and  time. 
The  slope  of  this  smoothing  operation  was  made  compatible  by  a  linear 
Interpolation  of  data  points  about  the  connecting  data  point.  The 
Infinite  space  solution  was  found  In  Erin gen  and  Suhubl  [20].  As  can  be 
seen  In  Figures  4.5  and  4.6  this  procedure  worked  well  and  Is  only 
necessary  at  or  near  the  free  surface. 


1-0  2j0  XO 

dimensionless  time.  C^/R 


Fig.  4.3  Green's  Function  and  Comparison  with  Payton 
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Fig.  4.6  Green's  Function  G22  anc*  Comparison  with  Eringen 
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V  ANTI-PLANE  STRAIN  WAVE  MOTION 

5.1  Introduction 

This  section  Is  specifically  concerned  with  developing  a  BEM 
technique  for  anti-plane  strain  using  a  time-dependent  approach  similar 
to  that  of  Cole,  et.al.  [15],  The  Green's  function  for  the  half-space 
geometry  will  be  selected  for  use  and  this  choice  will  eliminate  the  need 
to  discretize  the  free  surface  of  the  half-space  Itself,  thus  reducing 
still  further  the  number  of  required  unknown  variables.  Attention  will 
be  directed  to  SH-waves  propagating  In  a  half-space  with  cavities  or 
Inclusions.  The  wave  forms,  cavity  shape,  size  and  number  are  all 
arbitrary. 

5.2  Problem  Formulation 

The  direct  boundary  Integral  equation  method  will  be  used  to 
determine  the  solution  of  the  anti-plane  strain  case  for  transient 
elastodynamlc  waves  In  an  Isotropic,  homogeneous  half-space.  Recall  that 
anti -plane  strain  Is  characterized  by  a  displacement  field  of  the  form  ji 
*  {0,0,w(x-j  ,x2)}  .  The  non-zero  component  w  satisfies  the  classical 

wave  equation  from  (2.19),  l.e. 


w 


.11 


1 

77 


w. 


(5.1) 


where  c £  is  the  shear  wave  velocity,  and  the  body  forces  are  neglected. 

From  equation  (2.32)  the  Integral  representation  for  anti-plane 
strain  was  given  by 


00 

[GT(r0»to)"w(i:0’to)^-]dS(-o)dto 
o  as 


cw(r,t)  = 


(5.2) 


Equation  (4.5)  provides  the  proper  Green's  function,  l.e. 


H(t-t0-R,/c2) 


[(t 


-tJ2-R0Z/c2Z]'yi  +  [(t-t„)2-R,2/c,2]'/Z 


)  (5.3) 


1  '2 


and  the  associated  stress  field  Is  given  by  (4.6)  to  be 


1  * *  T  3" 


1  .  R0H(t-VVc2>  3Ro.  RlH(t-VRl/c2)  3R!  ,  ,C 
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Results  (5.3)  and  (5.4)  can  now  be  used  with  equation  (5.2)  to 
construct  a  boundary  Integral  equation  to  discretize  and  solve  for  the 
displacements  and  tractions  on  the  boundaries.  With  the  boundary  data 
determined,  equation  (5.2)  can  be  reapplied  to  find  the  displacements  and 
stresses  at  any  interior  points.  To  solve  for  the  stress  at  Interior 
points,  relation  (5.2)  Is  substituted  Into  (2.18)  yielding 


0  3S 


(5.5) 


5.3  Numerical  Discretization 

The  Integral  equation  (5.2)  must  now  be  discretized  In  some  manner  In 
order  to  model  a  physical  problem.  For  this  case,  the  boundary  to  be 
Integrated  over  was  divided  Into  elements  with  a  node  at  the  center  of 
each.  Each  Individual  element  Is  assumed  to  possess  a  constant 
displacement  and  a  constant  traction  with  respect  to  space.  Integration 
over  time  was  constant  for  the  displacement  kernel  and  linear  for  the 
stress  kernel.  Integration  over  a  particular  time  step  and  space 
Increment  requires  special  attention.  Rice  [45]  gives  a  description  of 
the  Integration  of  the  kernels  with  special  emphasis  on  the 


singularities.  This  was  essentially  the  same  scheme  used  by  Cole  [14]. 
It  should  be  pointed  out  that  it  is  possible  to  interpolate  the  variation 
of  the  dynamic  variables  over  an  element  with  any  order  of  shape  function 
(e.g.  constant,  linear,  quadratic,  etc.). 

For  a  receiver  point  on  the  boundary,  equation  (5.2)  thus  becomes 


m  H 


w  <r0*>3  dr0<tt0 

m  H  J  "  J 

=  I  E  /  /  G  (r  »t  )t?  dr  dt 

„_i  •  •,  +  r.  rJ  -o  0  J  0  0 

n-i  t  r j 


(5.6) 


where  H  is  the  total  number  of  elements,  and  M  Is  the  total  number  of 
time  Increments.  The  receiver  time  Is  incremented  by  m,  i.e.  t^  *  m^t, 
m  *  1,2,...M,  while  the  source  time  Is  Incremented  by  n,  tn  *  nA  t,  n  * 
1.2....M.  The  source  boundary  has  H  elements,  and  there  are  H  receiver 
elements.  Note  that  since  a  source  time  tn  cannot  be  greater  than  the 
receiver  time,  n  is  summed  only  up  to  m.  <j>j  and  ip ”  are  the 
Interpolation  functions  for  displacement  and  traction  respectively. 
Defining  the  discrete  kernels  or  Influence  coefficients 

IG™  -  /  /  Gdr  dt 

ij  t  r.  00 
n 

IKfjn  *{  '  ^r^^-ndr0dt0, 

n  j 

equation  (5.6)  can  be  written  as 


m  H  __  —  m  H  mn  n 

1  wT  +  E  E  IK1^  w"  =  E  E  IG'J'J 

2  1  n*l  j-1  J  n=l  j=l  1J  J 


(5.7) 


Equation  (5.7)  will  yield  a  set  of  H  linear  algebraic  equations  to  be 
solved  at  each  time  step.  The  final  result  of  this  will  be  the  tractions 
and  displacements  of  all  the  boundary  elements  at  all  the  time  steps. 


Although  there  are  2H^m2  possible  kernels,  not  all  need  to  be 
calculated  as  the  Green's  function  has  the  time  translation  property, 

G(r»t;r0,t0)  *  G(r,t+At;rQ,t0+At)  (5.8) 

and  the  causality  property 

G(r,t;r0,t0)  =  0  for  cg( t-tQ)< | r-rQ |  (5.9) 

These  properties  greatly  reduce  the  number  of  non-zero  kernels  in  a 
given  computational  problem.  Furthermore  it  is  interesting  to  note  that 
this  discretization  scheme  can  be  implicit  or  explicit  depending  on  the 
time  spacing  and  the  size  of  the  elements.  If  c2At£Ah/2,  where  h  is 
the  smallest  element  length,  the  coefficient  matrix  of  the  linear 
algebraic  system  to  be  solved  will  be  diagonal  and  the  discretization 
scheme  will  be  explicit.  If  c2At>  Ah/2,  the  matrix  will  no  longer  be 
diagonal  and  the  discretization  scheme  will  be  Implicit.  The  integral 
equation  Is  exact  so  any  errors  are  a  result  of  numerical  Integration. 
The  surface  of  the  body  Is  discretized  Into  small  line  segments  and  as 
one  would  expect  the  smaller  the  line  segments  the  more  accurate  the 
surface  representation  andthe  more  accurate  a  solution.  Integration  of 
the  Green's  function  In  time  has  dominant  errors  near  the  singularities. 
These  singularities  occur  when  the  direct  wave  or  reflected  wave  just 
reaches  the  receiver  point,  as  can  be  seen  In  the  equation  for  the 
Green's  function  (5.3).  As  discussed  in  Cole  [14],  smaller  time  steps  in 
general  will  give  better  accuracy.  The  problem  with  using  smaller  time 
steps  Is  that  one  will  need  a  larger  number  of  time  Increments  to  get  a 
suitable  displacement  output  and  there  is  a  buildup  of  roundoff  error 
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with  each  new  time  step. 

In  the  case  of  an  Interior  point,  the  displacement  and  stress  can  be 
solved  for  directly,  once  the  boundary  data  Is  known,  using  the  equations 

. .  m  H  mn  n  m  H 

w  =  Z  Z  IK™?  W?  +  Z  Z  IG1"?  t? 
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where,  p  Is  related  to  the  Interior  point  and  L  Is  the  normal  direction 


to  the  surface  the  stress  Is  on. 


5.4  Results 


An  example  problem  which  will  employ  the  solution  method  is  shown  In 
Figure  5.1,  and  contains  a  wave  source  and  a  stress  free  circular  cavity 
embedded  In  an  Isotropic,  homogeneous  half-space.  The  geometry  of  the 
source  was  a  four-sided  square  with  an  element  on  each  side  and  the 
cavity  was  discretized  Into  twelve  elements.  The  distance  between  the 
source  and  cavity  was  held  constant  (  i  -  0.1  m)  while  the  depth  of  the 
source  and  cavity  was  varied  (h  »  0.05,  0.075,  0.1  m).  The  size  of  the 
source  was  0.001  m  wide  and  the  cavity  was  0.01  m  in  diameter.  The 
matelal  properties  used  corresponded  to  a  shear  wave  speed  of  3,200  m/s 
and  a  shear  modules  of  8.1  x  10^  N/m^.  The  time  increment  was  At* 
2.5  x  10“®  s.  Each  configuration  was  run  with  the  same  Input  condition 
of  a  step  displacement  pulse  of  w0*  0.01  m  over  the  whole  source 
boundary  with  a  duration  of  one  time  Increment. 


(-f/2,h/2) 


(-e/2,h) 


(O.h/2) 


(SOURCE) 


Fig.  5.1  Schematic  of  Buried  Hole  Problem 


The  resulting  displacements  from  the  above  example  were  computed  at 
four  field  received  points  (shown  In  Figure  5.1).  Figure  5.2  Illustrates 
the  displacement- time  response  of  the  receiver  point  on  the  surface  of 
the  half-space  at  the  coordinates  (0,0).  There  Is  a  distinct  peak  as  the 
wave  reaches  this  point.  The  closer  the  source  to  the  receiver  the 
greater  the  magnitude  of  the  peak. 

In  Figure  5.3  the  receiver  point  Is  half  way  between  the  cavity  and 
the  free  surface  at  coordinates  (0,h/2).  In  addition  to  the  peaks  from 
the  direct  wave  there  are  now  peaks  from  the  reflection  from  the  free 
surface,  reflection  from  the  cavity  and  reflection  from  the  free  surface 
via  the  cavity.  Note,  at  h/z  -  0.5  the  events  of  the  reflection  from  the 
free  surface  and  reflection  from  the  cavity  should  occur  at  approximately 
the  same  time  according  to  ray  theory.  In  general  the  longer  the  ray 
path  the  smaller  the  peak  magnitude  of  the  event.  This,  however,  may  not 
be  the  case  near  the  free  surface  of  the  cavity  or  near  the  surface  of 


the  half-space.  Figure  5.4  corresponds  to  the  receiver  point  at 
coordlntes  (-t/2,h/2)  and  again  Illustrates  several  wave  arrivals. 

In  Figure  5.5  the  receiver  point  Is  located  between  the  source  and 
the  cavity  at  coordinates  (-z/2,h).  The  direct  waves  for  the  various  h/z 
cases  all  occur  at  the  same  time  with  essentially  Identical  magnitudes. 
For  the  h/ z  9  1.0  and  0.75  cases,  the  reflections  from  the  cavity  both 
occur  at  the  same  time  with  approxlmatey  the  same  magnitude.  Reflections 
from  the  cavity  for  the  h/z  *  0.5  case  are  associated  with  a  larger 
displacement  magnitude  because  of  the  presence  of  the  free  surface 
reflected  wave. 

These  computer  runs  were  made  on  a  Prime  750  minicomputer,  and  a 
typical  run  used  a  total  of  16  elements  with  4  output  points  and  40  time 
steps.  This  resulted  In  a  total  run  time  of  6.5  cpu(s). 
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Fig.  5.2  Results  of  Buried  Hole  Problem  at  (0,0) 
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Fig.  5.4  Results  of  Buried  Hole  Problem  at  (-£/2,h/2) 
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DIMENSIONLESS  TIME,  t/At 

Fig,  5.5  Results  of  Buried  Hole  Problem  at  (-H/2,h) 
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VI  PLANE  STRAIN  WAVE  MOTION 

|  6.1  Introduction 

This  section  is  concerned  with  developing  a  BEM  technique  for  plane 
strain  elastodynamics  using  a  time-dependent  Green's  function.  As  with 
|  the  previous  SH-wave  case,  the  Green's  function  for  the  half-space 

geometry  will  be  selected  for  use.  This  choice  will  again  eliminate  the 
need  to  discretize  the  free  surface  of  the  half-space  Itself,  thus 
j  reducing  the  number  of  required  unknown  variables.  Interest  will  be 

directed  to  SV-waves  and  P-waves  propagating  in  a  half-space  with 
cavities  or  inclusions.  The  wave  forms,  cavity  shape,  size  and  number 
are  all  arbitray. 


6.2  Problem  Formulation 

The  direct  boundary  Integral  equation  method  will  be  used  to 
determine  the  solution  of  the  plane  strain  case  for  transient 
elastodynamlc  waves  In  an  Isotropic,  homogeneous  half-space.  Plane 
strain  Is  characterized  by  a  displacement  field  of  the  form  jj  = 
{ufx-j.xg)  ,  v(x-|,X2),0}  .  The  equation  of  motion  from  (2.12) 

with  zero  body  forces  Is  given  by 

yV2u  +  (X+y)V(7.u)  =  pu  .  (6.1 ) 


L 


L 


From  equation  (2.30)  the  integral  representation  for  plane  strain  is 

00 

a  j  j  [Sl(r,t)*u(r,t)K*n]dS(r)dt  . 
o  3S 


(6.2) 


The  Green's  function  for  the  plane  strain  half-space  domain  was 


developed  In  section  4.3,  and  Is  given  by 

G,j  =  2 k  tH(t-r”/c1)U^H(t-tsp)UGJP  +  H(t-tps)U^ 

+  H(t-r"/c2)U^  +  H(t-tSps)H(r"/c2-t)U^S  +  H(t-r'/c1  )U^ 


where  corresponds  to  the  Incident  field,  and  is  given  by  Hudson 

1 J  13 

[24].  The  stress  field  associated  with  this  Green's  function  Is  from 
(4.20) 

Kijk  *  CW1J  Srs  +  M(i1r  6JS  *  61s  V)]Grk,s  •  ,6-4) 

Results  (6.3)  and  (6.4)  can  now  be  used  with  equation  (6.2)  to 
construct  a  boundary  Integral  equation  to  solve  for  the  displacements  and 
tractions  on  the  boundaries.  With  the  boundary  data  determined,  equation 
(6.2)  can  be  reapplied  to  find  the  displacements  at  any  interior  points. 

6.3  Numerical  Discretization 

The  Integral  equation  (6.2)  must  now  be  discretized  in  order  to  model 
a  physical  problem.  Similar  to  the  SH-wave  case,  the  boundary  to  be 
integrated  over  was  divided  into  elements  with  a  node  at  the  center  of 
each.  Each  Individual  element  Is  assumed  to  possess  a  constant 
displacement  and  a  constant  traction  with  respect  to  space  and  time. 
This  differs  slightly  from  the  anti-plane  strain  case  where  the 
displacement  was  constant  In  space  and  linear  in  time.  Because  the 
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half-space  Green's  function  Is  a  system  of  waves  (direct,  reflective,  and 
surface)  it  possesses  multiple  singularities  in  its  values  for  a  typical 
source/ receiver  pair  at  different  times.  Integration  over  a  particular 
time  increment  In  which  a  singularity  appears,  requires  special 
attention.  A  Gauss  four  point  quadrature  was  employed  to  integrate  both 
In  space  and  time.  This  allows  one  to  sample  points  over  the  integral 
other  than  the  endpoints  which  may  be  singular.  Rice  [45]  gives  a 
description  of  the  integration  of  the  kernals  with  special  emphasis  on 
the  singularities. 

For  a  receiver  point  on  the  boundary,  equation  (6.2)  becomes 
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n  J 


where  H,  M,  m,  n,  4,”  and  4)j  are  the  same  as  previously  defined  in 
section  5.3 

Defining  the  discrete  kernels  or  Influence  coefficients 

IGij  'idrodto 

rn  j 


IK?1?  *  /  /  K  dr  dt 

J  r  =  0  0 

n  m 


equation  (6.5)  can  be  written  as 
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Equation  (6.6)  will  yield  a  set  of  H  linear  algebraic  equations  to  be 
solved  at  each  time  step.  The  final  result  of  this  will  be  the  tractions 
and  displacements  of  all  the  boundary  elements  at  all  the  time  steps. 
Note,  this  closely  resembles  the  scheme  used  for  anti-plane  strain;  the 
difference  being  that  G  and  J(  are  now  second  order  and  third  order 
tensors,  respectively.  Also,  the  Integration  is  now  over  the  receiver 
points.  In  the  anti-plane  strain  case  the  Integration  was  over  the 

source  points. 

Although  there  are  8hV  possible  kernels,  not  all  need  to  be 
calculated  as  the  Green's  function  has  the  time  translation  property, 
(5.8)  and  the  causality  property  (5.9).  These  properties  greatly  reduce 
the  number  of  non-zero  kernels  In  a  given  computational  problem. 

Furthermore  It  is  again  interesting  to  note  that  this  discretization 
scheme  can  be  Implicit  or  explicit  depending  on  the  time  spacing  and  the 
size  of  the  elements.  If  Ci  At  <  h/2,  where  h  Is  the  smallest  element 

length,  the  discretization  scheme  will  be  explicit.  If  c-jAt>  Ah/2,  the 

discretization  scheme  will  be  Implicit.  This  method  works  better  for  an 
Implicit  scheme  where  At  c2  Is  approximately  equal  to  the  diameter  of 
the  largest  cavity.  When  At  c2  Is  less  than  the  diameter  of  the 
largest  cavity,  the  solution  blows  up  or  gets  larger  with  Increasing 
time.  There  appears  to  be  a  need  for  each  element  to  receive  Influence 
from  each  of  the  other  elements  on  a  continuous  surface. 

In  the  case  of  an  Interior  point,  the  displacement  can  be  solved  for 
directly,  once  the  boundary  data  Is  known,  l.e. 
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where,  p  Is  related  to  the  Interior  point. 

6.4  Resul ts 

In  section  4.4,  a  check  of  the  plane  strain  half-space  Green's 
function  was  made  by  comparing  the  result  with  the  Payton  [41]  and 
Erlngen  [20]  solutions.  As  a  che:k  of  the  total  BEM  scheme,  a  comparison 
was  made  between  Garvin's  analytical  solution  [22]  and  a  BEM  computer 
model.  Garvin's  work  provides  an  exact  solution  for  the  surface  response 
due  to  a  burled  dllatatlonal  line  source  of  Impulsive  time  dependence. 

In  trying  to  compare  these  two  solutions  It  should  be  pointed  out 
that  there  are  basic  unreasol vabl e  differences  between  the  Garvin 
solution  and  the  BEM  computer  model.  Garvin's  solution  Is  due  to  a 
burled  line  source  excitation  as  represented  In  Figure  6.1.  The  radial 


Fig.  6.1  Garvin  Problem  Configuration 


displacement  of  the  source  Is  given  by  Garvin  to  be 


at  H(t-r/c^ ) 
ZrttMrV)]1'2  ' 


(6.8) 


where  'a'  Is  a  constant  representing  the  source  strength  and  r  Is  the 
distance  between  the  source  and  receiver  point.  This  Is  taken  to  be  an 
approximate  representation  of  a  pulse  emitted  by  an  explosion,  and 
represents  a  sudden  singular  jump  followed  by  a  gradual  Incomplete 
recover,  see  Figure  6.2. 


Fig.  6.2  Garvin  Source  Displacement  Function 

In  contrast  the  BEM  model  Incorporates  a  burled  cavity  of  finite  size 
as  a  source  as  shown  In  Figure  6.3.  This  burled  source  was  modeled  as  an 
eight-sided  cavity  with  a  diameter  of  0.005  m  located  at  a  depth  of  0.1  m 


Fig.  6.3  BEM  Source  Configuration 


in  the  half-space.  On  this  burled  cavity  a  radial  traction  boundary 
condition  was  used  to  represent  the  source  excitation  In  Garvin's 
solution.  The  form  of  this  traction  function  was 


r6xl0 
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(6.9) 


and  Is  plotted  In  Figure  6.4.  Note,  the  function  Is  not  smooth  as  It  is 
discretized  to  be  constant  over  each  time  Increment. 

Therefore  In  order  to  alleviate  the  differences  between  the  actual 
Garvin  problem  and  the  BEM  model: 

1.  The  cavity  was  made  small  with  respect  to  its  depth. 
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2.  The  Garvin  solution  was  scaled  by  adjusting  the  source  strength 
parameter  'a1,  to  match  the  computer  solution  at  a  point  below 
the  surface  near  the  cavity.  A  comparison  was  then  made  at 
various  other  points  along  the  free  surface. 


Fig.  6.4  BEM  Boundary  Traction  Profile 


Figure  6.5  shows  the  vertical  displacement  solution  below  the  free 
surface  at  (0,7.5)  cm.  The  value  of  'a'  In  Garvin's  solution  was 
adjusted  to  be  1.3  x  10"^°cm^  to  match  the  BEM  solution.  The 
material  properties  correspond  to  a  P-wave  speed  of  5800  m/s,  an  S-wave 
speed  of  3350  m/s,  and  a  density  of  2500  kg/m^.  Figure  6.5  shows  both 
Garvin's  source  solution  and  the  BEM  solution  with  this  value  of  'a'  for 
comparison.  As  shown  in  Figure  6.6,  on  the  free  surface,  there  is  very 
good  comparison  directly  above  the  burled  source  at  (0,0).  Reasonably 
good  comparisons  also  exist  at  points  farther  away  from  the  source  as 
Illustrated  In  Figures  6.7-6.10. 


Fig.  6.5  Source  Comparisons  of  the  Vertical 
Displacement  at  (0,7.5)cm. 


Another  example  problem  which  will  employ  the  solution  method  is 
shown  In  Figure  6.11.  The  problem  (similar  to  the  SH-wave  example) 
contains  a  wave  source  and  a  stress  free  circular  cavity  embedded  in  an 
isotropic,  homogeneous  half-space.  The  geometry  of  the  source  was  a  four 
sided  square  and  the  cavity  was  modeled  as  a  twelve  sided  hole.  The 
distance  between  the  source  and  cavity  was  0.1  m  and  the  depth  of  the 
source  and  cavity  was  0.1  m.  The  size  of  the  source  was  0.001  m  wide  and 
the  cavity  was  0.01  m  in  diameter.  The  material  properties  used 
corresponded  to  a  P-wave  speed  of  5940  m/s  and  a  S-wave  speed  of  3,200 
m/s  and  a  density  of  7800  kg/m^.  The  time  increment  was  At=3.2  x 
10“®s.  This  configuration  was  run  with  the  input  condition  of  a  step 
normal  displacement  pulse  of  a*  0.001  m  over  the  whole  source  boundary 
with  a  duration  of  one  time  Increment. 


(0.0) 


Fig.  6.11  Schematic  of  Buried  Hole  Problem 


The  resulting  displacements  from  this  example  were  computed  at  four 
field  receiver  points  as  shown  in  Figure  6.11.  Figures  6.12-6.15  plot 
dimensionless  displacement  vs.  dimensionless  time.  The  dimensionless 
displacements  u/AQ  and  v/AQ  are  the  horizontal  and  vertical  displacements 
repsectively,  normalized  by  the  initial  step  displacement  magnitude.  The 
dimensionless  time  t/  At  Is  the  actual  time  divided  by  the  time 
Increment.  The  notation  in  these  four  figures  represent  the  events  that 
occur  as  particular  waves  pass  the  receiver  point.  These  waves  are 
classified  as 
Direct  Waves 

-  P:  (Dilatation  Wave) 

-  S:  (Shear  Wave) 

Free  Surface  Reflected  Waves 

-  PP:  (Incident  P-wave  and  reflected  P-wave) 

-  PS:  (Incident  P-wave  and  reflected  S-wave) 

-  SS:  (Incident  S-wave  and  reflected  S-wave) 

-  SP:  (Incident  S-wave  and  reflected  P-wave) 

Surface  Wave 

-  S-SP:  (Incident  S-wave  which  travels  along  the  surface  as 

P-wave) 

Cavity  Reflected  Waves 

-  CR-P:  (P-wave  reflected  from  the  cavity) 

-  CR-S:  (S-wave  reflected  from  the  cavity) 

Multiple  Reflected  Waves 


-SCR-P:  (P-wave  reflected  from  the  free  surface  and  then  from 
the  cavity) 


Figure  6.12  Illustrates  the  displacement-time  response  of  the 
receiver  point  on  the  surface  of  the  half-space  at  the  coordinates 
(0,0).  There  are  relatively  few  events  in  this  Figure  since  no  reflected 
waves  interact  with  this  point.  There  are  distinct  peaks  as  the  wave 
reaches  this  location  and  the  S- wave  and  surface  P-wave  act  together  at 
the  same  time.  Note,  because  the  ray  angle  makes  equal  angles  with  the 

X1  and  X£  axis,  the  horizontal  and  vertical  displacements  are 
approximately  the  same  from  the  time  the  P-wave  passes  until  another 
event  takes  place. 

In  Figure  6.13  the  receiver  point  is  half  way  between  the  cavity  and 
the  free  surface  at  coordinates  (0,5)  cm.  There  are  more  events  as 
reflections  from  the  free  surface  and  the  burled  cavity  are  observed. 
Note,  the  horizontal  displacements  are  relatively  larger  than  the 
vertical  displacements,  due  to  the  shallow  angle  between  the  source  and 
receiver  point.  Figure  6.14  corresponds  to  the  receiver  point  at 
coordinates  (-5,5)  cm.  and  again  Illustrates  several  wave  arrivals. 
Again  because  of  geometry,  the  horizontal  and  vertical  displacements  are 
approximately  equal  from  the  time  the  P-wave  passes  until  another  event 
takes  place.  In  Figure  6.15  the  receiver  Is  located  between  the  source 
and  the  cavity  at  coordinates  (-5,10)cm.  Note,  the  vertical  displacement 
Is  zero  until  reflections  occur. 

A  typical  computer  run  of  16  elements  with  4  output  points  and  30 
time  steps,  results  in  a  run  time  of  approximately  50  cpu(m). 
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VII  CONCLUSIONS  AND  RECOMMENDATIONS 


This  research  focused  on  the  feasibility  of  using  a  half-space, 
time-dependent  Green's  function  as  a  fundamental  solution  for  the 
construction  of  a  direct  BEM  code  to  solve  half-space  wave  propagation 
problems  for  In-plane  and  out-of-plane  motions.  With  this  half-space 
Green's  function,  the  need  to  discretize  the  free  surface  was  eliminat¬ 
ed.  The  time  dependent  approach  used  here  allows  the  causality  principle 
to  reduce  the  number  of  required  computations  and  permits  an  explicit 
time  stepping  scheme  to  be  constructed. 

The  present  state  of  the  method  and  code  has  some  limitations  and 
disadvantages.  Since  all  previous  time  steps  are  required  for  the 
solution,  large  amounts  of  computer  storage  are  necessary  to  use  the 
methods  developed.  Another  disadvantage  Is  that  while  the  method  is  good 
for  transient  problems,  application  to  the  steady-state  harmonic  problem 
reveals  solution  error  propagation.  This  Is  due  to  the  large  number  of 
time  steps  required  and  the  need  of  all  Information  from  the  previous 
time  steps.  Errors  In  the  method  arise  from  Integrating  the  Green's 
function  over  Its  singularity.  The  steady-state  harmonic  problems  have 
been  successfully  solved  using  BEM  with  Laplace  and  Fourier  transform 
methods. 

The  example  problem  containing  a  buried  source  and  a  stress  free 
circular  cavity  was  modeled  for  both  out-of-plane  and  In-plane  motions. 
For  this  problem  one  can  see  the  effect  of  a  much  more  complicated  wave 
phenomena  for  the  In-plane  (plane  strain)  case  when  compared  to  the 
SH-wave  case.  This  Is  due  to  the  effect  of  the  free  surface  as  It 
provides  mode  coversions  of  Incident  P-  or  SV-waves.  Because  of  these 


complexities,  the  time  to  compute  a  typical  run  was  considerably  longer 
for  the  plane  strain  problem. 

This  great  Increase  In  computation  time  over  the  SH-wave  case  is  due 
In  part  to  the  complex  nature  of  the  Green's  function.  To  solve  for  its 
values  for  a  typical  source/receiver  pair,  a  Newton  iterative  root 
finding  method  must  be  employed  (as  explained  In  section  4.3).  The  root 
finding  scheme  Is  also  used  to  find  the  ray  path  distances  as  described 
In  section  2.5.  There  Is  also  a  considerable  Increase  In  the  number  of 
components  of  the  Green's  function;  from  two  for  the  anti-plane  strain 
case  to  seven  for  the  plane  strain  problem.  In  addition  the  Integration 
Itself  Is  more  time  consuming  for  plane  strain.  In  the  anti-plane  strain 
case.  Integration  was  carried  out  using  a  differencing  scheme  from  Cole 
[14].  For  plane  strain  a  Gauss  four  point  quadrature  was  used  which 
means  the  Green's  function  Is  evaluated  16  times  for  each  space  increment 
and  time  step. 

In  general,  the  developed  methods  are  more  accurate  If  At  c2  Is 
approximately  equal  to  the  largest  cavity  size.  When  At  c2  is  less 
than  the  dimension  of  the  largest  cavity  the  solution  suffers  convergence 
problems  with  Increasing  time.  It  appears  that  each  element  should 
receive  Influence  from  each  of  the  other  elements  on  a  continuous 
surface.  Additional  numerical  experimentation  Is  needed  to  more 
accurately  assess  this  phenomena. 

The  results  obtained  here  for  these  plane  strain  and  anti -plane 
strain  problems  show  that  the  BEM  with  a  time-dependent,  half-space 
Green's  function  is  a  viable  method  for  elastodynamlc  half-space 
problems.  However  additional  research  Is  needed  In  order  to  develop  the 
technique  Into  highly  efficient  numerical  codes.  The  method  could  be 


extended  to  include  other  Interpolation  functions  (linear,  quadratic, 
etc.).  Although  Cole  [14],  has  Indicated  that  lower  order  functions 
would  be  best,  explicit  studies  of  this  point  should  be  made. 
Nonhomogeneous  continuum  problems  (layered  media)  could  be  investigated. 
Specifically  for  the  anti-plane  strain  case  one  could  use  more  Image 
sources  to  simulate  multiple  layers.  Anisotropic  or  nonlinear  continuum 
could  also  be  investigated,  using  the  general  method  developed,  providing 
the  half-space  Green's  function  could  be  found.  The  computability  of  the 
Green's  function  Is  very  important  on  the  overall  performance  of  the 
BEM.  The  current  state  of  the  code  developed,  could  certainly  be 
enhanced  by  further  improvements  of  the  plane  strain  Green's  function. 


-68- 


VIII.  CONTRACT  PUBLICATIONS 

1.  “Theoretical  Ground  Motions  Using  the  Boundary  Element  Method,"  Sadd, 
M.H.  and  Rice,  J.M.,  Proceedings  Fourth  Sensor  Technology  Symposium, 
Waterways  Experiment  Station,  Vicksburg,  MS,  1983,  (to  appear). 

2.  "Propagation  and  Scattering  of  SH-Waves  in  Semi-Infinite  Domains 
Using  a  Time  Dependent  Boundary  Element  Method",  Rice,  J.M.  and  Sadd, 
M.H.,  to  appear  In  Journal  of  Applied  Mechanics,  also  to  be  presented 
at  the  ASME  Winter  Annual  Meeting,  New  Orleans,  LA,  December  1984. 

3.  "El astodynamic  Boundary  Integral  Equation  Solutions  to  the  Half  Space 
Problem  Using  a  Time  Dependent  Green's  Function:  SH-Waves",  Rice, 
J.M.  and  Sadd,  M.H.,  Twentieth  Annual  Society  of  Engineering  Science 
Meeting,  University  of  Delaware,  1983. 

4.  "A  Note  on  Computing  El astodynamic  Full  Field  Displacements  Arising 
from  Subsurface  Singular  Sources",  Rice,  J.M.  and  Sadd,  M.H., 
Submitted  to  Mechanics  Research  Communications. 

5.  "In-Plane  Ground  Motions  Using  a  Time  Dependent  Boundary  Element 
Method",  submitted  to  International  Journal  of  Numerical  4  Analytical 
Methods  In  Geomechanics. 

6.  "Transient  Dynamic  Solutions  of  Some  Half  Space  Problems  by  the 
Boundary  Element  Method",  Ph.D.  Thesis,  J.M.  Rice,  University  of 
Rhode  Island,  May  1984. 


IX.  PARTICIPATING  SCIENTIFIC  PERSONNEL 


1.  Dr.  Martin  H.  Sadd,  Principal  Investigator. 

2.  John  M.  Rice,  Ph.D.  Graduate  Student,  Received  his  Ph.D.  In  May  1984 

3.  Krlshnaswamy  Kumar,  M.S.  Graduate  Student,  Employed  In  Summer  1982. 


BIBLIOGRAPHY 


1.  Abdul jabbar,  Z. ,  Datta,  S.K.,  and  Shah,  A.H.,  “Diffraction  of 
Horizontally  Polarized  Shear  Waves  by  Normal  Edge  Cracks  In  a  Plate," 
J.  Appl.  Phys.,  Vol.  54,  No.  2,  1983,  pp.  461-472. 

2.  Achenbach,  J.D.,  Wave  Propagation  In  Elastic  Solids,  North  Holland, 
Amsterdam,  1973. 

3.  Achenbach,  J.D.  and  Brlnd,  R.J.,  "Scattering  of  Surface  Waves  by  a 
Sub-Surface  Crack,"  J.  Sound  and  Vlb.,  Vol.  76,  1981,  pp.  43-56. 

4.  Achenbach,  J.D.,  Gautesen,  A.K.,  and  McMaken,  H.,  Ray  Methods  for 
Waves  In  Elastic  Solids,  Pitman,  London,  1982. 

5.  Alder,  B.,  Fernbach,  S.,  and  Rotenberg,  M.,  Methods  In  Computational 
Physics,  "Seismology:  Surface  Waves  and  Earth  0sc111at1ons7,r  VolY 
ll,  "Seismology:  Body  Waves  and  Sources,  Vol.  2,  Academic  Press, 
N.Y.,  1972. 

6.  Banerjee,  P.K.  and  Butterfield,  R. ,  Boundary  Element  Methods  In 
Engineering  Science,  McGraw-Hill,  England,  1981. 

7.  Berezin,  I.S.  and  Zhidkov,  N.P. ,  Computing  Methods,  Pergamon  Press, 
New  York,  1965. 

8.  Bouchon,  M.  and  K.  Akl,  "Discrete  Wave-Number  Representation  of 
Seismic-Source  Wave  Fields",  Bull.  Seism.  Soc.  Amer. ,  Vol.  67,  p. 
257,  1977. 

9.  Bouchon,  M. ,  "The  Importance  of  the  Surface  or  Interface  P-Wave  in 
Near-Earthquake  Studies",  Bull.  Seism.  Soc.  Amer.,  Vol.  68,  p.  1293, 
1978. 

10.  Bouchon,  M.,  "Discrete  Wave  Number  Representation  of  Elastic  Wave 
Fields  In  Three-Space  Dimensions",  J.  Geophys.  Res.,  Vol.  84,  p. 
3609,  1979. 

11.  Brebbla,  C.A.,  The  Boundary  Element  Method  for  Engineers,  Halsted 
Press,  New  York,  19781 

12.  Brebbla,  C.A.  and  Walker,  S. ,  Boundary  Element  Techniques  In 

Engineering,  Newnes-Butterworths,  London",  19801  “ 

13.  Cagnlard,  L.,  Reflection  and  Refraction  of  Progressive  Seismic  Waves, 
McGraw-Hill,  New  York,  1962. 

14.  Cole,  D.M.,  "A  Numerical  Boundary  Integral  Equation  Method  for 
Transient  Motion,"  Doctoral  Dissertation,  California  Institute  of 
Technology,  Pasadena,  Calif.,  1980. 

15.  Cole,  D.M. ,  Kosloff,  D.D.,  and  Minster,  J.B.,  "A  Numerical  Boundary 
Integral  Equation  Method  for  Elastodynamics  I,"  Bull.  Seism.  Soc. 
Amer.,  Vol.  68,  1978,  pp.  1331-1357. 


16.  Cruse,  T.A.,  “A  Direct  Formulation  and  Numerical  Solution  of  the 

General  Transient  El astodynami c  Problem  II,"  J.  Math.  Anal.  Appl., 
Vol.  22,  1968,  pp.  341-355. 

17.  Cruse,  T.A.  and  Rizzo,  F.J.,  "A  Direct  Formulation  and  Numerical 

Solution  of  the  General  Transient  Elastodynamic  Problem  I,",  J.  Math. 
Anal.  Appl.,  Vol.  22,  1968,  pp.  244-259. 

18.  Datta,  S.K.  and  El-Aklly,  N. ,  “Diffraction  of  Elastic  Waves  in  a 

Half-Space  1:  Integral  Representation  and  Matched  Asymptotic 

Expansions,"  Modern  Problems  In  Elastic  Waves,  Miklowitz,  J. ,  and 
Achenbach,  J.D.,  ed.,  Wiley,  New  York,  1978,  pp.  197-218. 

19.  de  Hoop,  A.T.,  "Representation  Theorems  for  the  Displacement  in  an 

Elastic  Solid  and  their  Application  to  Elasto-dynamic  Diffraction 

Theory,"  Doctoral  Dissertation,  Tech.  Hogesehole,  Delft,  1958. 

20.  Erinqen.  A.C.  and  Suhubl,  E.S.,  Elastodynamics,  Academic  Press,  New 
York,  1975. 

21.  Ewing,  W.M.,  Jardetsky,  W.S.,  and  Press,  F. ,  Elastic  Waves  in  Layered 
Media,  McGraw  Hill,  1957. 

22.  Garvin,  W.W.,  "Exact  Transient  Solution  of  the  Buried  Line  Source 

Problem,"  Proc.  Roy.  Soc.  London,  Vol.  234A,  1956,  pp.  528-541. 

23.  Graff,  K.F. ,  Wave  Motion  in  Elastic  Solids,  Ohio  State  Univ.  Press, 
1974. 

24.  Hudson,  J.A.,  The  Excitation  and  Propagation  of  Elastic  Waves, 

Cambridge  University  Press,  Cambridge,  1980.  ' 

25.  Kelly,  K.R.,  Ward,  R.W.,  Trellel,  S.,  and  Alford,  R.M. ,  "Synthetic 
Seismograms:  A  Finite  Difference  Approach,"  Geophysics,  Vol.  41, 
Feb.  1976,  pp.  2-27. 

26.  Ketter,  R.L.  and  Sherwood,  P.P.,  Modern  Methods  of  Engineering 

Computation,  McGraw-Hill,  New  York,  196'9.  “ 

27.  Kobtyashl,  S.  and  Nlshlmura,  N.,  "Dynamic  Analysis  of  Underground 
Structures  by  the  Integral  Equation  Method,"  4th  International 
Conference  on  Numerical  Methods  In  Geomechanics,  Edmonton,  Canada, 
1982. 

28.  Kobayashl,  S.  and  Nlshlmura,  N. ,  "Green's  Tensors  for  Elastic 
Half-Spaces,"  Memoirs  of  the  Faculty  of  Engineering",  Vol.  42,  Part 
2,  April  1980,  pp.  228-241. 

29.  Kupradze,  V.D.,  "Dynamic  Problems  in  Elasticity,"  Vol.  Ill, 
North-Holi  and,  Amsterdam,  1963. 

30.  Lamb,  H.,  "On  the  Propagation  of  Tremors  Over  the  Surface  of  an 
Elastic  Solid,"  Phil.  Trans.  R.  Soc.  London,  A.,  Vol.  203,  p.  1,  1904. 


-72- 


31.  Manolls,  G.D.,  “A  Comparative  Study  on  Three  Boundary  Element  Method 
Approaches  to  Problems  In  Elastodynamics,"  Num.  Meth.  In  Engineering, 
Vol.  19,  1983,  pp.  73-91. 

32.  Mansur,  W.J.,  and  Brebbia,  C.A. ,  “Formulation  of  the  Boundary  Element 
Method  for  Transient  Problems  Governed  by  the  Scalar  Wave  Equation," 
Appl.  Math.  Modelling,  Vol.  6,  August  1982,  pp.  307-311. 

33.  Mansur,  W.J.  and  Brebbia,  C.A.,  “Numerical  Implementation  of  the 
Boundary  Element  Method  for  Two  Dimensional  Transient  Scalar  Wave 
Propagation  Problems,"  Appl.  Math.  Modelling,  Vol.  6,  August  1982, 
pp.  299-306. 

34.  Mendelsohn,  D.A. ,  Achenbach,  J.D.,  and  Kerr,  L.M. ,  "Scattering  of 
Elastic  Waves  by  a  Surface-Breaking  Crack,"  Wave  Motion,  Vol.  2, 
1980,  pp.  277-292. 

35.  Mendelson,  A.,  "Boundary  Integral  Methods  in  Elasticity  and 
Plasticity,"  NASA  Tech.  Note  T.N.D.-7418,  173,  36  pp.,  1973. 

36.  Mlkhlln,  S.G.,  An  Advanced  Course  of  Mathematical  Physics, 

North-Holi  and,  19701!  ~~ 

37.  Mlsljenovlc,  D.M.,  "Boundary  Element  Method  and  Wave  Equation,"  Appl. 
Math.  Modelling,  Vol.  6,  August  1982,  pp.  205-208. 

38.  Mlsljenovlc,  D.M.,  "The  Boundary  Element  Method  In  Shock  Waves 

Analysis,"  Appl.  Math.  Modelling,  Vol.  6,  August  1982,  pp.  312-315. 

39.  Nlwa,  Y.,  Kobayshl,  S.,  and  Azuma,  N.,  "An  Analysis  of  Transient 
Stresses  Produced  Around  Cavities  of  Arbitrary  Shape  During  the 
Passage  of  Travelling  Waves,"  Memoirs  Fac.  Engng.,  Kyoto  Unlv.,  Vol. 
37,  Part  2,  1975,  pp.  28-46. 

40.  Nlwa,  Y.,  Kobayshl,  S.,  and  Fukui,  T.,  "Applications  of  Integral 
Equation  Methods  to  Some  Geomechanical  Problem,"  Numerical  Methods  In 
Geomechanics,  Desal,  C.S.,  ed.,  ASCE,  New  York,  1976,  pp.  120-131. 

41.  Payton,  R.G. ,  "Epicenter  Motion  of  an  Elastic  Half-Space  Due  to 

Burled  Stationary  and  Moving  Line  Sources,"  Int.  J.  Solids 

Structures,  Vol.  4,  1968,  pp.  28/ -300. 

42.  Reynolds,  A.C.,  "Boundary  Conditions  for  the  Numerical  Solution  of 
Wave  Propagation  Problems,"  Geophysics,  Vol.  43,  1978,  pp.  1099-1110. 

43.  Rice,  J.M.  and  Sadd,  M.H.,  "A  Note  on  Computing  Elastodynamic  Full 
Field  Displacements  Arising  from  Subsurface  Singular  Sources", 
submitted  to  Mechanics  Research  Comnunlcatlons. 

44.  Rice,  J.M.  and  Sadd,  M.H.,  "Propagation  and  Scattering  of  SH-Waves  In 

Semi- Infinite  Domains  Using  a  Time  Dependent  Boundary  Element 

Method,"  Journal  of  Applied  Mechanics,  to  appear. 


-73- 


45.  Rice,  J.M.,  “Transient  Dynamic  Solutions  of  Some  Half  Space  Problems 
by  the  Boundary  Element  Method",  Ph.D.  Thesis,  University  of  Rhode 
Island,  1984. 

46.  Sanchez-Sesma,  F.J.,  Herrera,  I.,  and  Aviles,  J.,  "A  Boundary  Method 
for  Elastic  Wave  Diffraction:  Application  to  Scattering  of  SH-Waves 
by  Surface  Irregularities"  Bull.  Seism.  Soc.  of  Amer. ,  Yol .  72,  pp. 
473-490,  1982. 

47.  Shah,  A.H. ,  Wong,  K.C.,  and  Oatta,  S.K.,  "Diffraction  of  Plane 

SH-Waves  in  a  Half-Space,"  Earthg.  Engg.  Str.  Dyn. ,  Yol.  10,  1982, 

pp.  519. 

48.  Shaw,  R.P.,  "Boundary  Integral  Equation  Methods  Applied  to  Wave 

Problems,"  Developments  in  Boundary  Element  Methods-I,  Applied 
Science  Pub,  London,  1979. 

49.  Sokol nikoff,  I.S.,  Mathematical  Theory  of  Elasticity,  Second  Edition, 
McGraw-Hill,  New  York,  1 956^ 

50.  Stone,  S.F.,  Ghosh,  M.L.,  and  Mai,  A.K.,  "Diffraction  of  Antiplane 

Shear  Waves  by  an  Edge  Crack,"  J.  Appl.  Mech.,  Yol.  47,  1980,  pp. 

359-362. 

51.  Telles,  J.C.F.  and  Brebbla,  C.A.,  "Boundary  Element  Solution  for 
Half-Plane  Problems,"  Int.  J.  Solids  Str.,  Yol.  7,  1981,  pp. 
1149-1158. 

52.  Wheeler,  L.T.  and  Sternberg,  £.,  "Some  Theorems  in  Classical 

Elastodynamics,"  Arch.  Rat.  Mech,  Anal.,  Yol.  31,  1968,  pp.  51-90. 


’•V 


APPENDIX 


COMPUTER  PROGRAMS:  BASIC  STRUCTURE 

Anti-Plane  Strain  Case 

Figure  A.l  Illustrates  the  flow  chart  for  the  anti-plane  strain 
program.  The  following  are  descriptions  of  the  contents  of  the  main 
program  and  subroutines. 

MAIN  PROGRAM:  Calls  three  subroutines  IMPUT.FMAT,  and  INTER. 

IMPUT:  Calls  and  defines  all  the  variables  necessary  to  solve  the 
problem  from  a  data  file.  These  variables  are  the  number  and  coordinates 
of  the  boundary  elements  on  all  cavities,  number  of  cavities,  the 
boundary  conditions  on  each  element,  and  the  number  and  coordinates  of 
the  Internal  points  at  which  a  solution  Is  desired.  Also  defined  are  the 
shear  wave  speed,  shear  modulus,  the  time  Increment,  and  number  of  time 
steps. 

FMAT:  Calls  subroutine  GF  which  will  integrate  the  individual 

kernels  for  displacement  and  traction.  FMAT  also  calls  subroutine  KOD 
which  will  provide  the  boundary  conditions  at  each  time  step.  A  system 
of  simultaneous  linear  algebraic  equations  are  then  formed  from  these 
Integrated  kernels  and  known  boundary  conditions.  The  unknown  boundary 
conditions  are  next  solved  for  by  the  simultaneous  equation  solver 
subroutine  SLNPD  using  Gauss  elimination.  Finally  subroutine  STORE  is 
called  which  stores  the  displacement  and  traction  history  for  all  surface 
elements. 


HAIM  PROGRAM 


Anti -Plane  Strain  Program  Structure 


6F:  Integrates  the  kernels  using  subroutines  ARBIT,  FGSH,  and  FKSH 
which  are  all  given  In  detail  In  Cole  [14], 

INTER:  Calculates  the  displacement  at  any  point  In  the  field  except 
cavity  surfaces. 

PLOTT:  Plots  the  displacement  vs.  time  for  the  field  points 

discussed  In  subroutine  INTER. 

Plane  Strain  Case 

Figure  A. 2  shows  the  flow  chart  for  the  plane  strain  program.  The 
following  are  descriptions  of  the  contents  of  the  main  program  and 
subroutines. 

MAIN  PROGRAM:  Calls  three  subroutines  IMPUT,  FMAT,  and  INTER. 

IMPUT:  Calls  and  defines  all  the  variables  necessary  to  solve  the 
problem  from  a  data  file.  These  variables  are  the  number  and  coordinates 
of  the  boundary  elements  on  all  cavities,  number  of  cavities,  the 
boundary  conditions  on  each  element,  and  the  number  and  coordinates  of 
the  Internal  points  at  which  a  solution  Is  desired.  Also  defined  are  the 
mass  density,  shear  wave  speed,  dilatation  wave  speed,  the  time 
Increment,  and  number  of  time  steps. 

FMAT:  Calls  subroutine  PSGF  which  will  Integrate  the  Individual 
kernels  for  displacement  and  traction.  FMAT  also  calls  subroutine  KOD 
which  will  provide  the  boundary  conditions  at  each  time  step.  A  system 
of  simultaneous  linear  algebraic  equations  are  then  formed  from  these 
Integrated  kernels  and  known  boundary  conditions.  The  unknown  boundary 
conditions  are  next  solved  for  by  the  simultaneous  equation  solver 
subroutine  MATSOL  using  Gauss  elimination.  Finally  subroutine  STORE  Is 
called  which  stores  the  displacement  and  traction  history  for  all  surface 
elements. 


PSGF:  Integrates  the  kernels  containing  the  Green's  function  and 
associated  stresses  using  a  Gauss  four  point  quadrature.  If  any  Green's 
function  does  not  decay  properly  after  all  the  wave  fronts  have  passed 
then  subroutine  MODD  Is  called  as  a  smoothing  operation  (as  discussed  in 
section  4.4). 

HUD'S:  HUD0,HUD1 ,HUD2,HUD3,HUD4,  and  HUD5  are  subroutines  that  calcu¬ 
late  Individual  components  of  the  Green's  function. 

HUDO:  Solves  for  the  Infinite  space  component. 

HUD1 :  Solves  for  the  PP  component. 

HUD2:  Solves  for  the  SP  component. 

HUD3:  Solves  for  the  PS  component. 

HUD4:  Solves  for  the  SS  component. 

HUD5:  Solves  for  the  SPS  component. 

STR'S:  STR0,STR1,STR2,STR3,STR4,  and  STR5  are  subroutines  that 
calculate  Individual  components  of  the  stress  associated  with  the  Green's 
function. 

STRO:  Solves  for  the  Infinite  space  component. 

STR1 :  Solves  for  the  PP  component. 

STR2:  Solves  for  the  SP  component. 

STR3:  Solves  for  the  PS  component. 

STR4:  Solves  for  the  SS  component. 

STR5:  Solves  for  the  SPS  component. 

NRSP:  Newton's  method  is  employed  to  solve  for  the  incident  S-wave 
and  reflected  P-wave  portion  of  the  Green's  function.  This  Is 

accomplished  by  computing  the  roots  of  complex  non-linear  equations. 

NRPS:  Newton's  method  Is  employed  to  solve  for  the  Incident  P-wave 
and  reflected  S-wave  portion  of  the  Green's  function.  This  Is 

accomplished  by  computing  the  roots  of  complex  non-linear  equations. 


NEWTPS:  Newton's  method  Is  used  to  evaluate  the  Incident  angle  of  a 
P- wave  and  the  reflection  angle  of  an  S-wave,  by  obtaining  the  solution 
of  two  simultaneous  non-linear  equations. 

NEWTSP:  Newton's  method  is  used  to  evaluate  the  incident  angle  of  an 
S-wave  and  the  reflection  angle  of  a  P-wave,  by  obtaining  the  sollution 
of  two  simultaneous  non-linear  equations. 

INTER:  Calculates  the  displacement  at  any  point  in  the  field  except 
on  cavity  surfaces. 

PLOTT:  Plots  the  displacement  vs.  time  for  the  field  points 

discussed  in  subroutine  INTER. 

A  complete  listing  of  the  computer  programs  may  be  obtained  by 
writing  to  Dr.  Martin  H.  Sadd,  Department  of  Mechanical  Engineering  and 
Applied  Mechanics,  University  of  Rhode  Island,  Kingston,  RI,  02881. 
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